---
title: "Figure 2"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## Note

The code in the files "ecolRxC_accuracy_IPF.html", "ecolRxC_accuracy_Thomsen_logit_no_aprox.html", "ecolRxC_accuracy_Thomsen_logit_Yule_aprox.html",  "ecolRxC_accuracy_Thomsen_probit_no_aprox.html" and "ecolRxC_accuracy_Thomsen_probit_Yule_aprox.html" must be previously run for this code to work. The package **tidyverse** is needed for this code to run.

## Introduction

This document contains the code to replicate Figure 2 in the paper.

## Data
```{r, eval = F}
# Data for Figure 2

output <- data.frame(specification = NULL, scale = NULL, measure = NULL, value = NULL)

# IPF
datos <- read.csv("summary_ecolRxC_solutions_IPF.csv", sep =";", stringsAsFactors = F)

EI.averages.logit <- apply(datos[, c(3:4,7:8)], 2, mean)
specification <- c("VTR-local", "VTR", "VTR-local-Yule", "VTR-Yule")
temp <- data.frame(specification, scale = rep("scale: logit", 4), 
                   measure = rep("EI: error index", 4), value = EI.averages.logit)
output <- rbind(output, temp)

EI.averages.probit <- apply(datos[, c(5:6,9:10)], 2, mean)
specification <- c("VTR-local", "VTR", "VTR-local-Yule", "VTR-Yule")
temp <- data.frame(specification, scale = rep("scale: probit", 4), 
                   measure = rep("EI: error index", 4), value = EI.averages.probit)
output <- rbind(output, temp)

EPW.averages.logit <- apply(datos[, c(11:12,15:16)], 2, mean)
specification <- c("VTR-local", "VTR", "VTR-local-Yule", "VTR-Yule")
temp <- data.frame(specification, scale = rep("scale: logit", 4), 
                   measure = rep("EPW: discrepancy index", 4), value = EPW.averages.logit)
output <- rbind(output, temp)

EPW.averages.probit <- apply(datos[, c(13:14,17:18)], 2, mean)
specification <- c("VTR-local", "VTR", "VTR-local-Yule", "VTR-Yule")
temp <- data.frame(specification, scale = rep("scale: probit", 4), 
                   measure = rep("EPW: discrepancy index", 4), value = EPW.averages.probit)
output <- rbind(output, temp)

EQ.averages.logit <- apply(datos[, c(19:20,23:24)], 2, mean)
specification <- c("VTR-local", "VTR", "VTR-local-Yule", "VTR-Yule")
temp <- data.frame(specification, scale = rep("scale: logit", 4), 
                   measure = rep("EQ: quadratic index", 4), value = EQ.averages.logit)
output <- rbind(output, temp)

EQ.averages.probit <- apply(datos[, c(21:22,25:26)], 2, mean)
specification <- c("VTR-local", "VTR", "VTR-local-Yule", "VTR-Yule")
temp <- data.frame(specification, scale = rep("scale: probit", 4), 
                   measure = rep("EQ: quadratic index", 4), value = EQ.averages.probit)
output <- rbind(output, temp)

# Thomsen_logit_FALSE
datos <- read.csv("summary_ecolRxC_solutions_Thomsen_logit_FALSE.csv", sep =";",
                  stringsAsFactors = F)

temp <- data.frame(specification = "ecolRxC", scale = "scale: logit", 
                   measure = "EI: error index", value = mean(datos$EI.ecolRxC.AVCR))
output <- rbind(output, temp)
temp <- data.frame(specification = "ecol-biN", scale = "scale: logit", 
                   measure = "EI: error index", value = mean(datos$EI.ecolRxC.Average))
output <- rbind(output, temp)


temp <- data.frame(specification = "ecolRxC", scale = "scale: logit", 
                   measure = "EPW: discrepancy index", value = mean(datos$EPW.ecolRxC.AVCR))
output <- rbind(output, temp)
temp <- data.frame(specification = "ecol-biN", scale = "scale: logit", 
                   measure = "EPW: discrepancy index", value = mean(datos$EPW.ecolRxC.Average))
output <- rbind(output, temp)


temp <- data.frame(specification = "ecolRxC", scale = "scale: logit", 
                   measure = "EQ: quadratic index", value = mean(datos$EQ.ecolRxC.AVCR))
output <- rbind(output, temp)
temp <- data.frame(specification = "ecol-biN", scale = "scale: logit", 
                   measure = "EQ: quadratic index", value = mean(datos$EQ.ecolRxC.Average))
output <- rbind(output, temp)


# Thomsen_logit_TRUE
datos <- read.csv("summary_ecolRxC_solutions_Thomsen_logit_TRUE.csv", sep =";",
                  stringsAsFactors = F)

temp <- data.frame(specification = "ecolRxC-Yule", scale = "scale: logit", 
                   measure = "EI: error index", value = mean(datos$EI.ecolRxC.AVCR))
output <- rbind(output, temp)
temp <- data.frame(specification = "ecol", scale = "scale: logit", 
                   measure = "EI: error index", value = mean(datos$EI.ecolRxC.Average))
output <- rbind(output, temp)


temp <- data.frame(specification = "ecolRxC-Yule", scale = "scale: logit", 
                   measure = "EPW: discrepancy index", value = mean(datos$EPW.ecolRxC.AVCR))
output <- rbind(output, temp)
temp <- data.frame(specification = "ecol", scale = "scale: logit", 
                   measure = "EPW: discrepancy index", value = mean(datos$EPW.ecolRxC.Average))
output <- rbind(output, temp)


temp <- data.frame(specification = "ecolRxC-Yule", scale = "scale: logit", 
                   measure = "EQ: quadratic index", value = mean(datos$EQ.ecolRxC.AVCR))
output <- rbind(output, temp)
temp <- data.frame(specification = "ecol", scale = "scale: logit", 
                   measure = "EQ: quadratic index", value = mean(datos$EQ.ecolRxC.Average))
output <- rbind(output, temp)


# Thomsen_probit_FALSE
datos <- read.csv("summary_ecolRxC_solutions_Thomsen_probit_FALSE.csv", sep =";",
                  stringsAsFactors = F)

temp <- data.frame(specification = "ecolRxC", scale = "scale: probit", 
                   measure = "EI: error index", value = mean(datos$EI.ecolRxC.AVCR))
output <- rbind(output, temp)
temp <- data.frame(specification = "ecol-biN", scale = "scale: probit", 
                   measure = "EI: error index", value = mean(datos$EI.ecolRxC.Average))
output <- rbind(output, temp)


temp <- data.frame(specification = "ecolRxC", scale = "scale: probit", 
                   measure = "EPW: discrepancy index", value = mean(datos$EPW.ecolRxC.AVCR))
output <- rbind(output, temp)
temp <- data.frame(specification = "ecol-biN", scale = "scale: probit", 
                   measure = "EPW: discrepancy index", value = mean(datos$EPW.ecolRxC.Average))
output <- rbind(output, temp)


temp <- data.frame(specification = "ecolRxC", scale = "scale: probit", 
                   measure = "EQ: quadratic index", value = mean(datos$EQ.ecolRxC.AVCR))
output <- rbind(output, temp)
temp <- data.frame(specification = "ecol-biN", scale = "scale: probit", 
                   measure = "EQ: quadratic index", value = mean(datos$EQ.ecolRxC.Average))
output <- rbind(output, temp)


# Thomsen_probit_TRUE
datos <- read.csv("summary_ecolRxC_solutions_Thomsen_probit_TRUE.csv", sep =";",
                  stringsAsFactors = F)

temp <- data.frame(specification = "ecolRxC-Yule", scale = "scale: probit", 
                   measure = "EI: error index", value = mean(datos$EI.ecolRxC.AVCR))
output <- rbind(output, temp)
temp <- data.frame(specification = "ecol", scale = "scale: probit", 
                   measure = "EI: error index", value = mean(datos$EI.ecolRxC.Average))
output <- rbind(output, temp)


temp <- data.frame(specification = "ecolRxC-Yule", scale = "scale: probit", 
                   measure = "EPW: discrepancy index", value = mean(datos$EPW.ecolRxC.AVCR))
output <- rbind(output, temp)
temp <- data.frame(specification = "ecol", scale = "scale: probit", 
                   measure = "EPW: discrepancy index", value = mean(datos$EPW.ecolRxC.Average))
output <- rbind(output, temp)


temp <- data.frame(specification = "ecolRxC-Yule", scale = "scale: probit", 
                   measure = "EQ: quadratic index", value = mean(datos$EQ.ecolRxC.AVCR))
output <- rbind(output, temp)
temp <- data.frame(specification = "ecol", scale = "scale: probit", 
                   measure = "EQ: quadratic index", value = mean(datos$EQ.ecolRxC.Average))
output <- rbind(output, temp)

write.table(output, "Figure2_data.csv", sep = ";", row.names = FALSE)
```

## Figure 2
```{r}
# Figure 2
library(tidyverse)
datos <- read.csv("Figure2_data.csv", sep =";", stringsAsFactors = T)

datos$specification <- ordered(datos$specification, 
                               levels = c("VTR-Yule", "VTR-local-Yule",
                                                 "ecol", "ecol-biN",
                                                 "VTR", "VTR-local",
                                                 "ecolRxC-Yule", "ecolRxC"))

p <- ggplot(datos, aes(x = specification, y = value)) +
      geom_bar(stat = "identity", fill = "azure3") +
      facet_grid(measure ~ scale, scales = "free_y") +
      scale_y_continuous(breaks = seq(0, 10, by = 5)) +
      theme_classic() +
      ylab(" ") + xlab("Procedure") +
     theme(axis.text.x = element_text(angle = 90, vjust = 0.3, hjust = 0.9),
           axis.title.y = element_text(angle = 90))


tiff("Figura2.tiff", height = 6.8, 
     width = 8, units = 'in', compression = "lzw", res = 300)
    p
dev.off()

p
```

